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Simulation of Spin Models in Multicanonical Ensemble with Collective Updates 
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We propose a Monte Carlo method which performs a random walk in energy space using cluster- 
like collective updates. By imposing that bond probabilities depend continuously on the micro- 
canonical temperature, we obtain dynamic exponents close to their ideal random walk values. The 
method proves remarkably powerful when applied to models governed by long-range interactions, 
where it straightforwardly combines with the efficient Luijten-Blote cluster algorithm to yield a 
dramatic reduction in the computation load. 

PACS numbers: 05.10.Ln, 64.60.Cn, 75.10.Hk 



Of the many methods dedicated to the study of spin 
models, Monte Carlo (MC) methods have now gained 
a prominent role. As regards models exhibiting cither 
first-order transitions or disorder-induced rough free en- 
ergy landscapes, canonical MC simulations are known to 
suffer, however, from supercritical slowing down pj, an 
exponential growth with the lattice size of the tunneling 
time between free energy minima that leads to unreli- 
able statistics. An efficient approach aimed at beating 
this limitation is the simulation in generalized ensembles 
0, 13 , in particular its multicanonical flavor initially pro- 
posed by Berg Q, 0] and independently by Lee , recon- 
sidered in the framework of transition matrix dynamics 
@ , and recently revisited by Wang and Landau \J\ . The 
key-idea here is to artificially enhance rare events corre- 
sponding to local maxima in the free energy, by feeding 
the Markovian chain with an auxiliary distribution W(E) 
best approximating the inverse of the density of states. 
Indeed, this was shown to reduce tunneling times from an 
exponential to a power law r ~ L z of the lattice size Q . 
Still, simulations in the multicanonical ensemble based 
on local updates yield dynamic exponents z which are 
substantially higher than their ideal random walk value 
z ~ D. In the case of long-ranged (LR) models, there is 
an additional hurdle owing to the very presence of long- 
range interactions which makes the computation of the 
energy — an essential ingredient of multicanonical algo- 
rithms — a very time consuming operation, namely, of 
0(L 2 ) complexity 0. 

In this Letter, we propose a Monte Carlo method which 
successfully tackles these two issues by performing sim- 
ulations in the multicanonical ensemble using collective 
updates. As opposed to previous approaches aimed at 
combining cluster updates with a multicanonical algo- 
rithm 9, 10, 111, our method relies on a straightforward 
cluster-building mechanism which hinges on the micro- 
canonical temperature of the current configuration in or- 
der to determine appropriate bond probabilities. We test 
the efficiency of the method on the two-dimensional q- 
state Potts model with nearest-neighbor (NN) interac- 



tions (q = 7, 10) and on its one-dimensional LR coun- 
terpart with l/r 1+<T interactions (q = 3,6, 12), with pa- 
rameters chosen so that a first-order regime is exhib- 
ited. In both test cases, analyses of tunneling rates 
show a very substantial reduction in the dynamic expo- 
nents, from e.g., z = 1.35(3) to z = 1.05(1) for the LR 
model with q = 6 and a = 0.7, and from z — 2.60(4) 
to z = 1.82(2) for the two-dimensional NN model with 
q = 7. We further demonstrate that our formulation 
makes it exceptionally straightforward to incorporate two 
acceleration schemes dedicated to LR models 

mm, 

which cut down the algorithm complexity from 0(L 2 ) to 
0(Lln(L)). Chains containing up to 2 16 spins were sim- 
ulated in a few days, whereas challenging such huge sizes 
with local updates would have demanded several months 
of intensive computation. 

Algorithm In the multicanonical method, one wishes 
to sample a flat histogram of the energy over a given 
energy range. The weight w(E) of a state of energy 
E is thus set to the inverse of the density of states, 
or more specifically, to an estimate of it obtained, in 
our case, using the Wang-Landau method 0. Denot- 
ing the microcanonical entropy as S(E), one may write 
w(E) = e~ s ( E \ A local- update algorithm consists in up- 
dating a single spin and accepting the attempted move 
with a probability given by min [l, e s ( Ea >~ s ( Eb >] , where 
E a and E^ stand for the energy of the initial and fi- 
nal states, respectively. In order to allow this algo- 
rithm to embody a collective-update scheme, we first 
rewrite the multicanonical weight w(E) as e~^^ E <j)(E) 
where (3(E) = dS(E)/dE is the inverse microcanoni- 
cal temperature. The very presence of a term having 
the same form as the canonical Boltzman weight pro- 
vides the means to reexpress the multicanonical weight 
as a trace over the bonds of a Fortuin-Kasteleyn random 
cluster 0, thus paving the way for a collective-update 
scheme. Although our algorithm may be equally well 
applied to other spin models, e.g., models incorporating 
disorder or exhibiting a continuous symmetry, we now 
consider, for the sake of clarity, a generalized ferromag- 
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netic spin model with a 1 q symmetry, whose energy reads 
E =- Ei<j J{\* - J'I)<W 3 ■ Here J(\i - j\) > and the 
Oi variables can take on integer values between 1 and 
q. Invoking the Fortuin-Kasteleyn representation of this 
model, we can write the multicanonical weight w(E) as 
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where the sum runs over all lattice bonds, a bond is ac- 
tive (inactive) whenever bij = 1 (0), and pu^ji(E) = 
e &(E)J(\i-j\) _ i, a collective-update step consists of two 
stages, namely first building a set of clusters from the 
current configuration, and then updating all clusters at 
once with an acceptance probability which ensures that 
detailed balance is preserved. The cluster construction 
can be carried out in exactly the same way as in the 
original Swendsen-Wang algorithm |15| . Starting from 
a configuration at energy E a and an empty bond set, 
we consider each pair of identical spins {er^, cr., } in turn, 
and place a bond with probability 1 — e~^ Ea ^ J " % ~^ . 
We then identify clusters of connected spins, and draw 
a new spin value at random for each cluster. Observing 
that a given bond configuration at energy E has a weight 
<j)(E)Yl l>0 pi(E) B W , where B(l) is the number of bonds 
of length I = \i — j\, we therefore accept the attempted 
cluster flips with the following acceptance rate: 



Wfu p (a — > b) = min < 1 
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Here, E\, denotes the energy of the new spin configura- 
tion. Since this acceptance rate is nothing but that of 
a Metropolis algorithm, detailed balance is trivially sat- 
isfied. In particular, if we consider a canonical simula- 
tion at inverse temperature /3o, we have (3(E) = /3 and 
4>(E) = 1; hence the acceptance rate is equal to 1, and we 
are back to the original Swendsen-Wang algorithm. It is 
crucial to underline that it is the microcanonical temper- 
ature, or equivalently the lattice energy, which entirely 
governs the cluster construction; indeed, for a given lat- 
tice configuration at energy E, bonds are placed as if the 
model were simulated at its micro-canonical temperature 
using a Swendsen-Wang algorithm. As a result, cluster 
bond probabilities change continuously as the lattice con- 
figuration walks along the available energy range of the 
random walk, so that, e.g., small clusters are built in the 
upper energy range and conversely large clusters in the 
lower energy range. Obviously, this mechanism entails 
determining (3(E) to sufficient accuracy (any departure 
from the ideal line resulting in poorer performances). We 
chose to compute (3(E) from the estimated density of 
states using spline interpolations, yet other means are 
available, e.g., one may rely on transition matrices [lfij |: 
this approach, though being slower, gave smoother esti- 
mates already at the very beginning of the simulation. 
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FIG. 1: Mean acceptance rate (W fu p ) as a function of the 
energy per spin and inverse microcanonical temperature (3(E) 
for the long-ranged six-state Potts chain model with a — 0.7, 
and L = 1024 spins. E and Edo denote the energy of the 
histogram peaks corresponding to the ordered and disordered 
phase, respectively. The finite size transition temperature T c 
and the 100% line are shown for convenience. 



We now examine two optimization schemes suited for 
long-ranged models. According to Eq.^ determining the 
acceptance rate of a cluster flip demands that we com- 
pute the energy of the new (attempted) lattice configura- 
tion beforehand. For long-ranged models, this represents 
an 0(L 2 ) operation, yet updating the lattice configura- 
tion in a collective way allows us to cut this complexity 
down to an O(LhiL) one by relying on an FFT imple- 
mentation of the convolution theorem 0, 0|. Still, it 
is crucial to note that this reduction is absolutely in- 
tractable with single-spin update implementations owing 
to the very reason that a single spin is updated at a 
time. For long-ranged spin models, the cluster-building 
process represents another exceedingly time-consuming 
operation, since at each MC step approximately L 2 pairs 
of spin are considered in turn for bond activation. When 
interactions decay with distance, a significant amount of 
time during the cluster construction is further wasted be- 
cause an overwhelming number of bonds considered for 
activation have only a negligible probability to be ac- 
tivated. Our formulation of the multicanonical weight 
w(E) in terms of a Fortuin-Kasteleyn mapping makes it 
straightforward, however, to build clusters using the effi- 
cient Luijten-Blote method based on cumulative bond 
probabilities ^| whereby, instead of considering each 
spin in turn for addition to a given cluster, it is the 
index of the next spin to be added which is drawn at 
random. The efficiency of this algorithm does not de- 
pend on the number of interactions per spin, and leads 
to a CPU demand which scales roughly as L. Our imple- 
mentation differs with that of ^3 essentially in that the 
cumulative bond probabilities now depend on the lattice 
energy through the microcanonical temperature, which 
is obviously constant over the whole cluster-construction 
process [l7j . 
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Numerical results We first discuss performance issues 
related to mean acceptance rates and tunneling times, 
for both the NN and the LR models. Indeed, as opposed 
to the Swendsen-Wang cluster algorithm, the acceptance 
rate of our algorithm (Eq. ^| is not trivially equal to 
unity, yet is tightly related to the efficiency with which 
the Markovian chain generates roughly independent sam- 
ples. An approximate analytical expression of the accep- 
tance rate when the initial and the final energies E a and 
E differ only by a small amount e is given to first order 
in e by Wfu p = min (1, 1 + A(E a )e), with 



A(E a ) = (3'{E a ) 
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The average energy is related to the average number of 
bonds of length I by (E) = - £ i>0 J(0^f (B(l)), 
which shows that (A(E)) — 0. Assuming a gaussian dis- 
tribution for A(E), it is easy to show that the variance 
(A(£) 2 ) is proportional to (3'(E) 2 and a term varying 
smoothly with E, whence 1 - (Wfu p ) oc \(3'(E)\e. As 
shown in Fig. ^ our numerical tests carried out on the 
six-state LR Potts chain with a — 0.7 show that (3(E) 
varies smoothly between the energy peaks of the ordered 
and disordered phases, which ensures that (Wfu p ) re- 
mains close to 1. The variance of A(£') increases when- 
ever E lies outside the range of phase coexistence, and, 
as is clearly visible, leads to a reduction in the acceptance 
rate. Still and all, it is worth underlining that the energy 
range of interest in the analysis of first-order phase tran- 
sitions spans an interval which is only moderately larger 
than the one corresponding to phase coexistence. In this 
respect, a mean acceptance rate remaining well above 
90% inside this range of energy represents already an im- 
provement of a factor 3 with respect to the standard mul- 
ticanonical approach where usual acceptance rates hardly 
exceed 30% .8] , let alone the crucial fact that a whole lat- 
tice sweep is now carried out in a single step. 

With regards to performance measurements at first- 
order transitions, tunneling times have so far been con- 
sidered one of the most meaningful measurement param- 
eters 

They are defined as one half of the 
average number of Monte Carlo update sweeps needed 
for the walk to travel from one peak of the energy his- 
togram to the other - where peaks are defined with 
respect to the finite-size transition temperature - and 
turn out to represent a good indicator of the interval 
between roughly independent samples. Results for the 
LR model with q = 3 and 6 are shown in Fig. Dy- 
namic exponents z were determined from a fit to the 
power law r e ~ L z . The collective-update algorithm 
yields z = 0.89(1) and z = 1.11(1) for q = 3, a = 0.4 and 
a = 0.6, and z = 1.05(1) for q = 6, a = 0.7. This rep- 
resents a substantial reduction with respect to the local- 
update implementation where we obtained, respectively, 
z = 1.13(2), z = 1.48(2) and z = 1.35(3); the improve- 
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FIG. 2: Tunneling times for the long-ranged Potts chain. Tri- 
angles and squares refer to the local- and collective-update 
algorithm, respectively. 



ment is even higher when one consider prefactors. For 
the NN model, simulations were performed over lattice 
sizes ranging from 16 x 16 to 256 x 256, giving z = 1.82(2) 
for q = 7 and z = 2.23(1) for q = 10. These estimates 
are much closer to the ideal value z ~ 2 expected from a 
random walk argument than those obtained with a local- 
update algorithm, namely z = 2.60(4) and z = 2.87(4), 
respectively. Additionally, they compare extremely well 
with those obtained with the multibond method and 
with Rummukainen's hybrid- like two-step algorithm 
although these approaches and ours differ markedly in 
the way clusters are constructed. 

In order to check that our algorithm did not pro- 
duce systematic errors, we computed transition temper- 
atures and interface tensions between coexisting phases 
for the NN model, for which exact results exist (|2JJ,|2lJ 
and references in For q = 10, we obtained 

T C (L) = 0.70699(5), 0.70300(2), 0.70278(1), 0.70164(1), 
and 0.701328(4) for L = 16, 30, 32, 64 and 128, where 
T c was determined from the location of peaks of the spe- 
cific heat. Following standard FSS theory at first-order 
transitions, we collapsed C V (T)/L 2 vs. (T — T C )L 2 over 
the four highest lattice sizes and found an infinite size 
temperature T c (oo) = 0.70123(5) in perfect agreement 
with the exact value 0.7012315 . . . The same procedure 
applied to q = 7 and L — 32, 64, 128 and 256 yielded 
T c (oo) = 0.77306(1) which again matches perfectly the 
exact value 0.7730589 . . . We estimated the interface ten- 
sion E from the histogram of the energy reweighted at a 
temperature where energy peaks have the same height, 
namely, 2E = — L~ x ln(P m i„), where P m in denotes the 
minimum of the histogram between the two energy peaks, 
and the peak heights are normalized to unity. Our algo- 
rithm allowed us to determine E with a four-digit pre- 
cision for sizes up to L = 256 and nonetheless rather 
modest statistics (of order 10 7 sweeps). For the seven- 
state NN model, we obtained 2E = 0.0336(6), 0.0294(1), 
0.02631(8) and 0.02384(9) for L = 32, 64M28 and 256; 
a linear fit of the form E ~ E(oo) + c/L 22] performed 
over the three largest sizes (i.e., for L above the disor- 
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FIG. 3: CPU time per MC step and per spin for the long- 
ranged Potts chain. Triangles indicate typical CPU times for 
the local-update algorithm, irrespective of q and a. Filled 
squares refer to the collective-update algorithm, where for 
q — 3 estimates were determined by averaging over a — 0.4, 
0.5 and 0.6. 



dered phase correlation length £ ~ 48 |2l() yielded the 
infinite size value 0.02230(11), still above the exact value 
0.020792, yet closer to it than estimates reported in sev- 
eral previous studies 0, 0, 0] . We note in passing that 
this discrepancy may be very well attributed to the influ- 
ence of higher order terms in the vicinity of L ~ £, since 
retaining the two largest sizes only would yield a closer 
value of 0.02137(20). 

Finally, we discuss CPU demand performances in the 
case of LR models. Assuming a decently efficient algo- 
rithm implementation, this indicator gives a rough ac- 
count of the algorithm complexity. Figure |21 sketches 
averages of the CPU (user) time per MC step and per 
spin. Small fluctuations might be attributed to the ef- 
fect of CPU caches differing in size. While for the local- 
update implementation the demand in CPU per spin 
grows linearly with the number of spins, it is roughly 
constant over a fairly large range of lattice sizes with the 
collective-update algorithm. Moreover, the local-update 
implementation is outperformed already at sizes of sev- 
eral hundreds spins, with nonetheless an increased foot- 
print for higher q due to the correspondingly higher num- 
ber of FFT's to be computed. This clearly demonstrates 
the breakthrough that this new method brings about for 
the numerical study of LR models, drawing in particular 
highly precise tests of finite-size scaling within computa- 



tion range. 

Our method is easily generalized to other spin models 
for which a cluster representation is available. Amongst 
promising candidates are disordered models, which are 
known to exhibit rugged energy landscapes and high dy- 
namic exponents; related implementation details and re- 
sults are reported in a distinct work [l7j] . 
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